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Abstract. The traditional multiconfiguration Hartree-Fock (MCHF) and configura- 
tion interaction (CI) methods are based on a single orthonormal orbital basis. For 
atoms with many closed core shells, or complicated shell structures, a large orbital 
basis is needed to saturate the different electron correlation effects such as valence, 
core- valence correlation and correlation within the core shells. The large orbital ba- 
sis leads to massive configuration state function (CSF) expansions that are difficult 
to handle, even on large computer systems. We show that it is possible to relax the 
orthonormality restriction on the orbital basis and break down the originally very 
large calculations to a series of smaller calculations that can be run in parallel. Each 
calculation determines a partitioned correlation function (PCF) that accounts for a 
specific correlation effect. The PCFs are built on optimally localized orbital sets and 
are added to a zero-order multireference (MR) function to form a total wave function. 
The expansion coefficients of the PCFs are determined from a low dimensional gener- 
alized eigenvalue problem. The interaction and overlap matrices are computed using 
a biorthonormal transformation technique (Verdebout et ai, J. Phys. B: At. Mol. 
Phys. 43 (2010) 074017). The new method, called partitioned correlation function 
interaction (PCFI), converges rapidly with respect to the orbital basis and gives total 
energies that are lower than the ones from ordinary MCHF and CI calculations. The 
PCFI method is also very flexible when it comes to targeting different electron corre- 
lation effects. Focusing our attention on neutral lithium, we show that by dedicating 
a PCF to the single excitations from the core, spin- and orbital-polarization effects 
can be captured very efficiently, leading to highly improved convergence patterns for 
hyperfine parameters compared with MCHF calculations based on a single orthogonal 
radial orbital basis. 

By collecting separately optimized PCFs to correct the MR function, the variational 
degrees of freedom in the relative mixing coefficients of the CSFs building the PCFs are 
inhibited. The constraints on the mixing coefficients lead to small off-sets in computed 
properties such as hyperfine structure, isotope shift and transition rates, with respect to 
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the correct values. By (partially) deconstraining the mixing coefRcients one converges 
to the correct limits and keeps the tremendous advantage of improved convergence 
rates that comes from the use of several orbital sets. Reducing ultimately each PCF 
to a single CSF with its own orbital basis leads to a non-orthogonal configuration 
interaction approach. Various perspectives of the new method are given. 
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1. Introduction 

The electron correlation energy of an atom has been defined by Lowdin as 
the difference between the exact nonrelativistic energy eigenvalue of the electronic 
Schrodinger equation and the energy of the single configuration state function (CSF) 
approximation, commonly called the Hartree-Fock energy. In line with this definition 
we think of electron correlation effects as those manifesting themselves beyond the 
Hartree-Fock approximation. It is useful to subdivide further and consider both static 
(nondynamical) and dynamical correlation [2]. Static correlation arises from near- 
degeneracies of the Hartree-Fock occupied and unoccupied orbitals. Systems with 
significant static correlation are poorly described by a single CSF and are said to have 
a strong multireference character. Dynamical correlation is due to the reduction in the 
repulsion energy related to the electron-electron cusp and is a short-range effect. 

Accurate description of electron correlation remains a major challenge in atomic 
structure calculations. To meet this challenge a number of different methods have 
been developed such as many-body perturbation theory (MBPT) [3l H], combinations 
of configuration interaction and many-body perturbation (CI+MBPT) [5], [6], [7], and 
coupled cluster (CC) [HI El theories. Different kinds of variational methods have also 
been used, and one may specially note Hylleraas-type calculations, that explicitly include 
the interelectron distance r^ in the construction of the wave function [HI [121 1131 E] • 
In quantum chemistry, variational complete active space self-consistent field (CASSCF) 
methods are quite successful for describing small and medium-size molecules, but are 
not sufficient when dynamical correlation must be included [15]. The latter are treated 
through second-order perturbation theory using a single or multireference state as 
the zero-order approximation. Combined variational multireference and second order 
Moller-Plesset perturbation calculations have also been applied very successfully by 
Ishikawa and co-workers [161 HZ] to obtain accurate transition energies for a number of 
atomic systems. In this paper we will critically examine variational multiconfiguration 
methods, such as multiconfiguration Hartree-Fock (MCHF) combined with configuration 
interaction (CI). 

Multiconfiguration methods are quite general, and can be directly applied to excited 
and open-shell structures across the whole periodic table. By including the most 
important closely degenerate CSFs to form a multireference (MR) expansion, the static 
correlation is efficiently captured. Dynamical correlation is accounted for by adding, 
to the MR expansion, CSFs obtained by single (S) and double (D) excitations from 
the CSFs in the MR to an increasing set of active orbitals. The CSFs generated in 
this way build a space that we refer to as the correlation function (CF) space and it 
is convenient to think of multiconfiguration expansions as something built from CSFs 
in the MR space and in the CF space. Due to restrictions in the Racah or the Slater 
determinant algebra underlying the construction of the energy expressions, the orbitals 
are usually required to be orthonormal. Such an orthonormal orbital basis is not very 
efficient for larger systems. Let us consider an atomic system, for example Ca, with 
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several closed shells. To describe the dynamic correlation in the Is shell we tailor 
an orbital set for which some of the orbitals should have a large overlap with the Is 
radial orbital. Due to the orthogonality restrictions of the orbitals, the correlation in 
the 2s shell needs to be described in terms of the previous radial orbitals, tailored 
for describing correlation within the Is shell, as well as some new radial orbitals that 
are overlapping with the 2s orbital etc. To capture the dynamic correlation between 
electrons in all the different shells, the orbital basis needs to be extended to a large 
number of orbitals for each symmetry, leading to massive CSF expansions [18j. This 
is in effect a scaling wall [19] that has been difficult to get around. In practice the 
electrons in the atom are considered as either core or outer valence electrons resulting 
in valence- valence, core-valence, and core-core types of SD excitations with the latter 
often neglected. Another general problem with variational methods is that they are 
entirely based on the energy functional, and properties not strongly coupled to this 
functional may be inadequately described by the resulting wave function. As an example 
we consider the hyperfine interaction. The CSFs that are responsible for the important 
spin- and orbital-polarization effects are relatively unimportant for the total energy, and 
thus the orbital basis from the variational calculation may be spatially localized in such 
a way that the above effects are not captured. Alternatively, a very large orbital basis 
is needed to achieve convergence for these properties, leading to CSF expansions that 
grow unmanageably large. 

The present work is an extension of a previous study of correlation energy 
in beryllium [20]. Based on a fast transformation technique, originally proposed 
by Malmqvist and collaborators [211 122], we show that it is possible to relax the 
orthonormality restriction on the orbital basis and use several mutually non-orthogonal 
orbital basis sets that are better adapted to the short range nature of the dynamical 
correlation. The gained freedom also makes it possible to tailor an orbital basis 
for capturing effects weakly connected to energy, improving convergence properties 
of atomic properties other than the energy. Partitioning the CF space into several 
subspaces, and using different orbital sets optimized for the different partitions, may be 
one way around the scaling wall associated with single orthonormal orbital sets. 

2. Partitioning the MCHF problem 

Starting from the non-relativistic Hamiltonian for an iV-electron system 



the multiconfiguration Hartree-Fock (MCHF) approach determines an approximate 
wave function of the form 

M 




N r 
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in which a CSF, $,(7^ LS'"), belongs either to the MR space or to the CF space. All 
CSFs have a given parity (vr) and LS symmetry and they are built from a common basis 
of one-electron spin-orbitals 

(j){nlmims) = ^P{nl; r)Yi^X^,<^)xms, (3) 

where the radial functions P{nl] r) are to be determined [22] ■ For the approximate 
wave function ([2]), the integro-differential MCHF equations have the form 

Yinl ; r)] - - e„,.,} P(n/ ; r) 

= '^X{nl;r) + Y,^ni,n'iP{n'l;r) (4) 

for the unknown radial functions [23]. The equations are coupled to each other through 
the direct Y and exchange X potentials and the Lagrange multipliers tniyi- The 
Lagrange multipliers force the radial orbitals to be orthonormal within the same / 
subspace. Under these conditions the configuration state functions are orthonormal 

= 6,j. (5) 

The mixing coefficients appearing in the expansion over CSFs also enter in the form of 
the potentials and are determined by solving the configuration interaction (CI) problem 

He = Ec, (6) 

with Hij = being the Hamiltonian matrix and c = (ci, C2, . . . , CmY the column 

vector of mixing coefficients. For a given set of mixing coefficients, the equations (jl]) 
are solved by the self-consistent field (SCF) procedure. The SCF and CI problems 
are solved, one after the other, until convergence of both the radial functions and the 
selected Cl-eigenvector is achieved. 

The strong coupling between the CSF expansion and the resulting optimized orbital 
basis is well known [21]. In the variational multiconfiguration approach indeed, the 
orbitals adapt spatially to account for the specific correlation effect targeted by the 
tailored expansion. In the present work, we investigate the possibility of breaking 
down the computational task into subtasks by partitioning the CF space into different 
subspaces, each targeting a specific correlation effect such as valence- correlation, core- 
valence or correlation within shells in the core, and performing separate MCHF 
calculations for each expansion built on the MR space and a CF subspace. In a final 
step the wave function is expanded in a basis consisting of CSFs from the MR space and 
functions built in each of the CV subspaces. The expansion coefficients are obtained 
by computing the Hamiltonian and overlap matrices and solving the corresponding 
generalized eigenvalue problem. The computation of the matrix elements between 
functions in the different subspaces, and this is the crucial point, is made possible 
by the biorthogonal transformation [2T| [22]. 

The above scheme offers various advantages, resembling the "Divide and Conquer" 
strategy - i) from the computational point of view, smaller subtasks can be run in 
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parallel, - ii) the resulting orbital basis sets are better adapted for capturing efficiently 
electron correlation, with the hope of getting a satisfactory accuracy for the desired 
property before reaching the scaling wall, - iii) the coupling between the subspaces 
reduces to a reasonably small dimension eigenvalue problem. 

3. The PCFI approach 

The efficiency of the method with a partitioned CF space was shown in our ffist 
paper [20], when targeting the total energy of the ground state of neutral beryllium. A 
multi-reference (MR) expansion 



limited to the major contributions to valence correlation, including the near-degenerate 
ls^2s^ and ls^2p^ configurations of the Layzer's complex, was corrected by three 
Pair Correlation Functions |Avv)) l^cv), |Acc) built by allowing single and double 
excitations from specific subshells of the MR configuration state functions to a given 
orbital active set (AS), and specifically tailored to describe the valence (VV), core- 
valence (CV) and core (CC) correlation effects. The final wave function 



yielded a lower energy than the traditional MCHF method based on a very large CSF 
expansion. 

In the present work, we generalize the approach by introducing the ^^Partitioned 
Correlation Functions" instead of ^^Pair Correlation Functions", preserving the PCF 
acronym, but allowing more flexible building rules for each PCF. For describing 
unambiguously the CSFs content of a given PCF, we first define a "pure" PCF as a 
CSFs expansion containing only one kind of excitations, ie. single (S), or double (D), or 
triple (T),. . .excitations. We write such a Partitioned Correlation Function A{j}, where 
the subscript {i} specifies the set of occupied shells that are excited to a given active set. 
For instance, a PCF including only double excitations from the ni/i and 712/2 subshell^j] 
is written Kn^i^n2i2- ^ PCF could also be "hybrid" if containing different kinds of 
excitations. Such a PCF may be seen as a superposition of "pure" PCFs and is written 
as A{j}_{j}_ , where for each family of excitations {z}, {j}, . . ., the letters sequence 
appearing in the subscript specify the n/-labels of the MR electrons that undergo the 
excitations. It is clear that a Partitioned Correlation Function in our approach does not 
fit with the usual definition of a pair correlation function [251 126] . 

According to this notation, a PCF representing single excitations from uili and 
double excitations from the uili and 71,2/2 subshells, is written as 



m 




(7) 



^(7 LS^)) = |^^^^^(7 LS^))+aYY |Avv) + acv |Acv) + acc |Acc)(8) 





X Note that nili and 712^2 may refer to equivalent or non-equivalent electrons. 
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where the first summation corresponds to all possible single excitations from the riili 
shell for each CSF belonging to the MR and the second one, to all possible double 
excitations from the ni/i and 7x2/2 shells for each CSF belonging to the MR. One can 
then rewrite the beryllium ground state wave function (|8]) as 

+q;vv \K-vv) + acv \Msv) + acc |Ais_isis) 

where v stands for any occupied valence subshell of the MR set. In order to avoid 
having the same CSF in two different partitioned correlation functions, we consider the 
partitioning of the configuration space into disjoint sets. Another possible partition 
satisfying this property could be 

1^(7 = |^^^(7 

+a'yy |A 

v—vv 

where the single excitations (Is — )■ nih) have arbitrarily moved from the CC to the 
CV correlation function subspace. In the most general case, the MR function is 
corrected by p PCFs 

|vl/(7 LSn) = |^"'^(7 LSn) + E |A,) , (10) 

1=1 

each of the PCFs corresponding to a given partition of the CF space: 

dim (A) 

|A) = E 4 i^t) ■ (11) 

We use the notation \E'^ for the function consisting of the MR function and one of the 
correcting PCFs A 

m dim(A) 

Lsn) = J24 i^^^'i^^ Lsn) + E 4 ■ (12) 

i=i j 
In our approach, this function is obtained by solving the corresponding MCHF 
equations (j4j)-(l6]) to optimize the A-PCF orbital set and mixing coefficients. Such a 
calculation that optimizes the MR eigenvector composition {af} with orbitals frozen 
to the MR-MCHF solution ([7]), the mixing coefficients {cf} and the A-PCF radial 
functions, will be referred as a MR-PCF calculation. As far as the notation is concerned, 
we will underline when necessary the orbitals that are kept frozen during the self- 
consistent-field process. Solving the MCHF problem ( IT2l) for each Aj [i = 1,2, ... ,p) 
produces p mutually non-orthonormal one-electron orbital sets. Each of the orbital sets 
will be optimally localized for the correlation effect described by the corresponding PCF 
expansion. Assuming the CSFs of the MR and CF spaces orthonormal and = 1, 

we have 

m dim(h) 

= (vi/MR|vi/MR^ + (A|A) = + Yl (4^' = 1 ' 
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revealing that (A|A) 7^ 1. To keep a natural interpretation of the PCF weights, we 
renormalize each PCF according to 

1 



A 



A . 



The PCFI approach consists in regrouping the m components of the MR space and 
the p CF subspaces in an a priori low-dimension interaction matrix to get a compact 
representation of the total wave function 

rn p 

i=i j 

where the mixing coefficients {oj} and {bj} are obtained by solving the generalized 
eigenvalue problem of dimension {M = m + p) 

He = ESc . (14) 

The corresponding Hamiltonian matrix may be explicitly written as 



H 



($MR|^|$MK^ 



MR\ 



($MR|^|$MK^ 



MR\ 



/<|.MR|^|$MR\ 

\ m I I m / 



($MR 
/$MR 

\ m 


l^lAi) 
l^lAi) 


(Ail^lAi) 



($MR|/j|Ap) 



^T\H\Ap) 



\ 



{Ai\H\Ap) 



\ 



{Ap\H\Ap) 



(15) 



The matrix dimension M is simply given by the sum of the number of CSFs belonging 
to the MR (m) and the number of PCF functions (p). The overlap matrix has the 
same structure, with a value of unity on the diagonal thanks to the renormalization 
(|Aj) — )■ |Aj)). It can be obtained by merely substituting the H operator, appearing 
in each matrix element, by the unit operator 1. It reduces to the unit matrix for 
specific PCF-building rules as shown in Appendix 1. All plain-line blocks involve 
orthonormal orbitals and the construction of the matrix elements between the CSFs 
in the blocks is based on fast angular integration methods developed by Gaigalas et 
al. [27t [28] . This holds not only for the diagonal blocks, but also for blocks coupling 
the CF and the MR spaces since we do not allow the MR orbitals to vary in the MR- 
PCF MCHF calculations of (^^. The building of all other blocks of the CF space, 
surrounded by a dashed line, involve non-orthogonal orbitals arising from independent 
MCHF calculations and requires therefore the use of biorthonormal transformations 
before the traditional methods for angular integration can be applied (see section 6 
of 120] )■ By solving this compact eigenvalue problem, we showed in [20] that accurate 
total energies can be obtained. 
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However, as we already mentioned in that work, some variational freedom in the 
coefficients is lost by the fact that solving f lT4|) does not allow relaxation in the relative 
weights within each PCF. The latter are indeed fixed linear combinations, 

|A) = 5^c^|$.) (16) 

k 

and as such, we will refer to |A) as a constrained CSFs expansion in the sense that 
the orbitals and the expansion coefficients are not allowed to change. The coefficients 
{c^} will be called the constrained coefficients. The effect associated with this loss of 
flexibility on a property, that we are investigating in the next section, will be qualified 
as the constraint effect. 

4. The constraint effect 

We were expecting that this constraint would have a minor impact on all expectation 
values but extending our previous work ^20j to ls^2s2p of Be brought a surprise for 
spectroscopic properties other than the total energy. For this state, we optimized the 
following MR (m = 6) function 

MR = ls^{2s2p, 2s3p, 2p3s, 2p3d, 3s3p, 3p3d ^P°} , 

that we corrected by three PCFs functions [p = 3) targeting respectively valence, core- 
valence and core correlation and built on single- and double-excitations (SD) from the 
multi-reference set 

\<il{ls^2s2p ip°)) = I^MR^ipo^^ 

+avv \^v-vv) + «cv l^isv) + «cc l^is-isis), 

where v stands for any orbital of the valence MR set, ie. v = {2s, 2p, 3s, 3p, 3d}. After 
solving the MR-PCF problem ( !T2|) for the three A's, we solved the (M = m + p = 9) 
eigenvalue PCFI problem. We will use the acronyms SD-MR-PCFI and SD-MR-MCHF, 
respectively, for the present Partitioned Correlation Function Interaction and for the 
conventional multi- reference MCHF calculations. Figure [H shows the striking difference 
in behavior between the resulting energy and specific mass shift (SMS) parameter. 
While the total energies obtained with the SD-MR-PCFI and the SD-MR-MCHF models 
converge to the same limit, the two approaches give rise to an unexpected off-set on the 
SMS parameter. 

For a deeper understanding of the ins and outs of the PCFI approach and fully 
appreciating its advantages, we moved to a smaller system: neutral lithium and its 
spectroscopic properties. For the lithium ground state, the Hartree-Fock approximation 
is rather good and the monoreference (m = 1) ls^2s "^S can be taken. We apply the 
PCFI method using two {p = 2) PCFs: -i) the ffist one targeting single and double 
excitations from the core (Is) orbital and denoted Ais-i^i^ , - ii) a second one, A2s-is2s, 
targeting single excitations from the 2s valence shell and double excitations from the core 
(Is) and valence (2s) orbitals. The size of the PCFI matrix is small (M = 3). Figure [2] 
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SD-MR-PCFI {ls,2s,2p,3s,3p,3d} 
SD-MR-DPCFI {ls,2s,2p,3s,3p,M} 
SD-MR-MCHF 




6 7 8 
")7iaT value 



10 11 



0.448 
0.446 
0.444 
0.442 
0.44 
0.438 
0.436 
0.434 
0.432 



SD-MR-PCFI {Is, 2s, 2p. 3s, 3p, 3d} 
SD-MR-DPCFl {Is, 2s, 2p, 3s, 3p, M} 
SD-MR-MCHF 
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Figure 1. Be second excited state. Comparison of the convergence trends between 
the SD-MR-PCFI, SD-MR-MCHF and SD-MR-DPCFI (see section E]) approaches for 
the total energy (left figure) and the SMS parameter (right figure). The horizontal 
lines correspond to the results of Komasa and Rychlewski [29] . 



illustrates the MCHF radial functions of the (n = 5) active set resulting from the two 
MR-PCF equations applied to Ais^i^i^ and A2s-is2s- As it was observed in [201 121] for 
other systems, one can realize from this figure that a given PCF orbitals set specifically 
favors the region of the space occupied in the reference by the electrons undergoing 
excitations. 

We compare the expectation values of other operators than the Hamiltonian, i.e. 
the specific mass shift and the hyperfine interaction parameters [30], evaluated by the 
two SD-(mono)-PCFI and SD-(mono)-MCHF methods in figure El The two curves illus- 
trate the impact of the "constraint effect" on three different properties: the total energy, 
the specific mass shift (S'sms) and the contact term (ctcont)- As mentioned in the end 
of section |3l the origin of this effect is the hindrance to free variation in the expansion 
coefficients and in the orbitals. Even if the pre-optimized PCF orbital sets are fixed in 
the PCFI diagonalization step, the associated orbital constraint effect is expected to be 
small in comparison to the constraint on the mixing coefficients. The lack of variation 
in the orbitals is indeed somewhat compensated by the use of separated PCF orbital 
sets that makes the number of radial functions for a given active set, larger in the PCFI 
approach than in the MCHF method. Moreover the number of correlation layers used 
for a given PCF is probably large enough to reach saturation. The strongest limitation 
is likely to come from the fact that for each PCF, the mixing coefficients appearing 
in ( [T6l) are kept frozen in the interaction step that leads to the final many-electron wave 
function expression ( !T3|l . By freezing these mixing coefficients, we inhibit the expression 
of any indirect effects. The relative weights of the CSFs within each PCF are indeed 
already fixed by the frozen coefficients {cf} obtained from each independent MR-PCF 
MCHF optimization of ( !T2|) . each one targeting a specific correlation component, and 
no possibility is offered to these coefficients to capture the higher-order PCF-coupling 
effects. 
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Figure 2. The two thick hnes correspond to the spectroscopic Hartree Fock Is (no 
node) and 2s (one node) orbitals of Li Is^ 2s '^S. The other lines represents the radial 
functions of the correlation orbitals of the n = 5 active set. The comparison between 
the top (core-valcnce MCHF) and bottom (core-core MCHF) figures illustrates the 
contraction of the first few correlation orbitals when going from h2s-is2s to Ais_isis 
(mono)-PCF calculations. 
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Figure 3. Convergence of the absolute total energy, the specific mass shift and the 
hyperfine contact parameters for the ground state of neutral lithium. The agreement 
found between the SD-MCHF and SD-PCFI curves for the total energy is not observed 
for the two other properties. We recover the consistency with the traditional method 
by deconstraining the wave function (SD-DPCFI). 



Engels [3T] studied the influence of various excitation classes on ab initio calculated 
isotropic hyperfine constants and showed how each class gives a direct contribution, 
but that there are also important secondary effects on the contributions from the 
other classes due to relative changes of mixing coefficients. The importance of these 
indirect effects on the CSF weights that affect computed properties quite substantially is 
confirmed in the present work. Back to the beryllium example, Tabled] that collects the 
radiative data for the resonance El transition of beryllium illustrates that the constraint 
effect can also be significant for properties involving more than one state. Comparing 
the gf values from the SD-MR-MCHF and SD-MR-PCFI calculations we see that the 
latter converges much more rapidly, but that there is an off-set of a little less than 1 %. 
Also, the agreement between the length and velocity forms is not as good for the PCFI 
method as for the ordinary MCHF. 



A PCFI approach for describing electron correlation in atoms 



13 



Table 1. Line and oscillator strength of the — resonant line for the beryllium 
atom. 
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1.37491 


9 






42580.53 


10.6363 


10.6301 


1.37570 


1.37490 


10 






42577.57 


10.6369 


10.6322 


1.37569 


1.37508 


CAS {fljYiaxlmax — 


= 9i) 




42588.71 


10.6234 


10.6333 


1.37431 


1.37559 


Chung and Zhu [! 


m 




42559.20 






1.374 




Fleming et al. [33 


I 




42593.44 






1.3745 


1.3759 


Komasa and Rychlewski [29] 


42560.52 










Glowacki and Mi^ 


;dalek 


m 


42670.69 






1.340 





Kramida and Martin [35] 
Irving et al. [36] 
Schnabel and Kock [37l 



42565.35 



1.40 ±0.04 
1.34 ±0.03 
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5. Deconstraining Partitioned Correlation Functions 

In the PCFI method, the expansion coefficients for the CSFs in the PCF are constrained 
(locked) so that there is no possibility of relative changes due to the interaction with 
other PCFs. To recover this variational freedom, the PCFs can be deconstrained by 
transferring hj CSFs from the j-th PCF to the basis and at the same time setting their 
weights to zero, i.e. extending the PCFI space to 

{$r, • • • , $r } U{^.} {^r, • • • , $r i u^^i' • • • > ^i'^'}- (i^) 

The superscript d for the PCF A"' indicates a renormalized c?e-constrained Partitioned 
Correlation Function whose weights of the transferred CSFs have been set to zero. The 
many-electron wave function expansion becomes 

m P ( f^j I 

i=i j=i [k=i J 

where the expansion coefficients are obtained from a higher dimension (M = m + 
X]j=i(^i + 1)) priori generalized eigenvalue problem. The size of each block in (ITSll 
involving at least one PCF is growing accordingly to the number of deconstrained CSFs. 
In the limit of the completely deconstrained case (ie. hj = dim{Aj) Vj ), we regain 
full variational freedom in the coefficients, with the advantage that each CSF brings its 
tailored orbital basis. M is then at his maximum value, ie. the total number of CSFs, 
and the wave function ( JTSll will be referred as being "deconstrained" . It is strictly equiv- 
alent to a configuration interaction problem in the CSF space built on various mutually 
non-orthonormal orbital sets. 



For solving the eigenvalue problem ( IT^ and for building efficiently the interaction 
matrices associated with the selected operators in the basis of CSFs and deconstrained 
PCFs spanning the wave function f lT8|) . we modify the original way of presenting the 
biorthonormal transformation [22] to evaluate Oij = (Aj|0| Aj), where O is the Hamilton 
or unit operator. Following the original formalism, we perform a biorthonormal 
transformation 

to express the original left and right hand side PCFs in the new CSF bases {$^1 and 

rii rii 

k=l k=l 
rij rij 

|A,) = E^'I*') = E^"^I'^')' (20) 
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where the counter-transformed eigenvectors {c/^} and {df} ensure the invariance of the 
total wave functions. Given the matrix representation O of an operator O 

= i^imi) , (21) 

the matrix element between these PCFs is written as 

(A^io|A,) = J2(^irdi{^imi) = (d^)*6d^- (22) 

k,l 

where d is the column vector of counter-transformed mixing coefficients. Note that each 
CSF expansion should be closed under de-excitation (CUD) for allowing the biothonor- 
mal transformation [22l [20] . By strictly following this methodology, we may think that 
we should apply one biorthonormal transformation for each matrix element associated 
with any off-diagonal sub-matrix block involving two non-orthogonal orbital sets. How- 
ever, as the counter-transformation process is fixed by the overlap between the original 
spin-orbital bases, i.e. {4>k\'PD — ^km^ show in the Appendix 1 that one can eval- 
uate the whole sub-matrix block by performing a single biorthonormal transformation 
treating simultaneouly the counter-transformation of all the elements constituting the 
block-basis. 

Using this strategy, relaxing the PCFI constraint to any degree becomes possi- 
ble. The price to pay is the increase of the size of the Partition Correlation Func- 
tion Interaction problem f lT^ . In the limit of the completely deconstrained case (ie. 
hj = dim{Aj) Vj) , the PCFI approach is strictly equivalent to a configuration in- 
teraction problem in the original complete CSF space. We will use the label DPCFI 
for this "deconstrained" approach. Since each CSF could be built, if worthwhile, on 
its own orbital basis, without any radial orthogonality constraints with the other ones, 
this DPCFI approach is equivalent to a general configuration interaction problem in 
non-orthogonal orbitals. 

The (D)PCFI procedure can be summarized as follows: 

(i) Perform a HF/MCHF calculation for the mono-/multi- reference wave function ([7]), 

(ii) Freeze the orbitals belonging to this MR space and perform p separate MR-PCF 
MCHF calculations ( !T2l) for the different Partitioned Correlation Functions, 

(iii) (Optional) Deconstrain each PCF by transferring the desired CSFs from the CF to 
the MR basis (see ([IT])), 

(iv) Build the Hamiltonian and other relevant operators interaction matrices ( |T5i) by 
performing the biorthonormal transformations, if necessary, using the weight matrix 
formalism (see Appendix 2), 

(v) Solve the eigenvalue problem ([T^ for getting the many-electron wave function ( [T3l) . 

(vi) Compute the desired property with the PCFI eigenvector (s). 

As it is shown in Appendix 1, step (ii) allows to replace the overlap matrix equal of ([T4ll 
by the unit matrix. 
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6. Applications of the DPCFI method 

6.1. Solving the constraint problem 

This "Deconstrained Partitioned Correlation Function Interaction" (DPCFI) strategy 
has been applied to both lithium and beryllium by solving the eigenvalue problem in 
the original CSF basis instead of the PCF one. Each CSF comes with the orbital basis 
associated with the PCF from which it comes from. In both cases, as illustrated by 
figure [Hand figure [3] for beryllium and lithium respectively, the DPCFI values converge 
to the MCHF results, recovering all the indirect effects and keeping the tremendous ad- 
vantage of a much faster convergence. This corroborates the fact that the problematic 
PCFI-MCHF discrepancy observed for the S'sms and contact hyperfine parameters is due 
to the constraint effect. On the other hand, total energies are only slightly improved in 
the (PCFI — DPCFI) transition. A close analysis of the results displayed in Figure [3] 
shows that for n = 10, the relative differences (in absolute value) between the PCFI and 
DPCFI results are respectively 1.65 10~^ %, 0.18 % and 1.33 % for the total energy, the 
specific mass shift and the contact term parameters. These numbers demonstrate that, 
compared with the total energy, the specific mass shift and the contact term are much 
more sensitive to the constraint effect. 

6.2. The beryllium resonance line 

For the beryllium resonance transition (see tabled]), passing from SD-MR-PCFI to SD- 
MR-DPCFI, the line strength is modified by about 0.7 % in the length formalism and 
1.1 % in the velocity formalism, bringing the two forms in much closer agreement. The 
SD-MR-DPCFI is now in excellent agreement with the most extensive MCHF calcula- 
tions available (n = 9 CAS). As the absolute energy of each level, the transition energy 
is much less sensitive to the constraint effect since it undergoes a modification of less 
than 0.02 %. Note that both the SD-MR-PCFI and the SD-MR-DPCFI methods pro- 
vide a transition energy at the n = 6 level which is comparable to the SD-MR-MCHF 
n = 10 result. As already pointed out in our first paper [20], the richness in terms of 
radial functions of the many-electron wave functions adopting different PCF definitely 
leads to a higher rate of convergence for all properties. 

It is probably worthwhile to review briefiy the status of the available oscillator 
strength values for this Be I resonance line reported in table [H Experimentally, the 
situation has evolved since the publication in 1996 of the theoretical results of Flem- 
ing et al. [33] of gfi = 1.3745 in the length form and gfy = 1.3759 in the velocity form. 
Irving et al. [36] revised the old beam-foil measurement of Martinson et al. |i38j, by in- 
cluding the cascade corrections through the ANDC analysis that increased the oscillator 
strength from gf = 1.34 ± 0.05 1.40 ± 0.04, in nice agreement with the theoretical 
predictions. An independent measurement has been realized by Schnabel and Kock [37] 
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using the cascade-free laser induced fluorescence method, yielding the original beam-foil 
result but with a smaller uncertainty, namely 1.34 ±0.03. Amazingly, this value escaped 
to the attention of the authors of a rather complete compilation [39]. Table [1] shows 
that the present SD-MR-DPCFI results confirm the previous theoretical value [33], with 
a difference between the two gauges that is much smaller than the accuracy indicators 
reported by Fuhr and Wiese [IQ]. The small overlap between the two most recent ex- 
perimental values [36| |37] call for further investigations on the experimental side. 

6.3. The CAS-DPCFI approach in lithium 

In this section we present results obtained for selected spectroscopic properties involving 
the two lowest states of neutral lithium, ie. ls^2s ^5 and ls^2p ^P°. The Li I ground 
state (D)PCFI calculations that were presented in section [Hand figure [3] to investigate 
the constraint effect are limited in the sense that triple excitations are systematically 
omitted in a SD-mono-reference calculation. In order to obtain more accurate many- 
electron wave functions, we adopted the complete active space (CAS) list of CSFs 
combined with the DPCFI approach that remains manageable for a three-electron 
system. For illustrating the great flexibility of the DPCFI approach, we investigate two 
different strategies - i) treating core excitations globally and ii) separating the single 
and double core-excitations and dedicating a PCF to core-polarization. 

6.3.1. A global core description Like in our first model, we use the Hartree-Fock (HF) 
solution as the zeroth-order wave function for the two lowest states of this three electron 
system. In this approach, we split the CF space made of single, double and triple 
excitations (SDT) in three different subspaces {p = 3) defining the following three 
PCFs: 

• one for taking care of the inner-shell correlation between the two Is electrons 

|Ai,_i,i,) = \ls^2l ^U) + \ls2l n'l' ^U) 

n'l' 

+ Y '^") , (23) 

n'l',n"l" 

• a second one associated to the inter-shell correlation between the Is and the 21 

{I = s or p) 

\K2i^is2i) = \ls^2l ^L^) + 5^ lls^ n'l' ^L^) 

n'l' 

+ Y \^sn'l'n"l" ^L^) , (24) 

n'l',n"l" 

• a third and last one including the pure triple excitations 

\Aisis2i) = \ls^2rL^)+ Y \n'l'n"l"n"'l"' ^L^) . (25) 

n'l',n"l",n"'l"' 
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The above notation for the three PCFs apphes to both ls^2s "^S or ls^2p states, 
with L = I. 



For the first two PCFs ( 123|) and fl24|) . we optimize all the correlation orbitals, 
freezing the Is and 21 (= 2s or 2p) to the HF solution of the mono-reference. This 
strategy is inadequate for Aisis2« since it only contains triple excitations that do 
not interact with the reference CSF. The optimization of the corresponding orbital 
set becomes then more tricky. In the previous paper [20j, we chose to use a SD- 
multireference to include triple excitations. In the present work, we dedicate a specific 
PCF to these. We first define an "extended" SD expansion for a reference set built on 
the [n = 2, 3) shells and optimize it by allowing variations in the correlation orbitals 
only. This expansion opens an indirect interaction between the triple excitations and the 
reference CSF. For capturing these higher-order effects, we optimize the n > 3 orbitals 
during the MR-PCF procedure. The n = 4 layer is therefore the first one that effectively 
represents three-electron excitations. Figure H] illustrates for the two states that the 
DPCFI convergence is faster than the traditional CAS-MCHF approach based on the 
same CSF expansions. For a given orbital active set, the corresponding total energy 
value is indeed systematically below the CAS-MCHF result. Since the angular content 
of the wave function (maximum /-value for in the one-electron basis) is identical for 
both methods, we conclude that the DPCFI method captures more efficiently electronic 
correlation for a given atomic system. A similar improvement is a priori expected for any 
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Figure 4. Convergence of the total energy with respect to the largest principal 
quantum number of the active sets for the ground state ls^2s of lithium on the left 
and for the first excited state ls^2p on the right. The reference values (dotted 
lines) correspond to the absolute energy values obtained by Yan et al. [41] . 

other spectroscopic property. Figure [5] presents the convergence pattern of the contact 
term for the ground state and the electric quadrupole parameter of the first excited state 
of neutral lithium. As it clearly appears, the hyperfine parameters are not converging 
as smoothly as the total energy. It is well known that the relevant expectation values 
are extremely sensitive to single excitations [21] and it is worthwhile to attempt another 
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Figure 5. Convergence of the hyperfine contact parameter of the ground state (on 
the left) and of the electric quadrupole hyperfine parameter of the first excited state 
(on the right) with respect to the largest principal quantum number of the active sets. 
The reference value for the contact Uc parameter (dotted line) corresponds to the value 
obtained by Yan et al. [32]. 



6.3.2. A Partitioned Correlation Function dedicated to core-polarization For describing 
more accurately the hyperfine interaction, we split the Ais„i5i5 PCF in two subspaces 

Ai,_isis ^ Ais + Aisis • (26) 

The Ais PCF function focusing on the single excitations is dedicated to capture core- 
polarization effects. The many-electron wave function is then written as the reference 
function corrected by four {p = 4) different PCFs with their own orbital set: the two 
first, associated to the core- valence and to the triple excitations ( l25l) . as described 
in the previous subsection, completed by 

• a third one that takes care of the single excitations of the Is shell 

|Ai,) = 115^2/ ^L") + 11^2/ n'l' ^U) , (27) 

n'V 

• a last one associated to the double excitations of the Is shell 

|Ai,i,) = \ls^2l ^L^) + \2l n'l'n'T ^L^) . (28) 

n'V,n"l" 

Excitations considered in (127|) describe spin-polarization, for both the "^S and 
states, since single excitations can break the singlet spin coupling between the 
two core electrons. For the first excited state, these excitations can also break the 
angular coupling associated with orbital-polarization. It is well known that the hyperfine 
parameters are sensitive to these excitations [13] and some improvement is expected in 
their evaluation thanks to the splitting The results are presented in figure [HI 

figure [7] and figure [HI By comparing figure H] and figure |6l it is obvious that the 
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decomposition (12^ does not affect the total energy value. The DPCFI method still 
captures correlation more efficiently than the traditional MCHF calculations. The 
interesting improvement appears for the different hyperfine parameters. Figure [7] and 
figure [8] illustrate their progressive convergence, respectively for the ground and the first 
excited states. The use of the orbital set tailored for capturing the spin- and orbital 
core-polarization (CP) enhanced beautifully the convergence pattern of all the hyperfine 
parameters. The resulting trends are much smoother than those of the global core 
approach and the ordinary MCHF (see figure [5l). All oscillations disappeared and we 
reach reasonably well- converged values around n = 5. Lithium is a small atom and it is 
possible to enlarge an ordinary orbital basis to get converged values for all quantities. For 
larger atoms with more complicated shell structures it is, to set things into perspective, 
often not possible to extend the radial orbital basis very much due to a rapidly growing 
number of CSFs and here the fast convergence of the CP-DPCFI method, together with 
the fact that orbital sets for different shells can be optimized independently of each 
other, represents a major improvement in the general methodology. 
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Figure 6. Convergence of the absolute energy with respect to the largest principal 
quantum number of the active sets for ls^2s ^S* (left) and ls^2p (right) in neutral 
lithium. The reference values (dotted line) correspond to the results obtained by Yan 
et al. gl]. 

The values of the corresponding total energies, S'sms and hyperfine parameters are 
reported in tables [2] and |31 As in the figures, we compare the DPCFI values with 
Yan et al. [HI |32] results using Hylleraas-type variational method. The remaining 
difference between both studies may be attributed to the slow angular convergence rate 
[El - Ei.i = 0{l + 1/2)-^) of the (l/ria) angular development [B]. 

Table H] presents the transitions energies, line strengths and weighted oscillator 
strengths obtained using the CAS-CP-DPCFI method. This table illustrates the slow 
convergence of the oscillator strength in both gauges. It leads to a remaining gap 
between the length and the velocity gauges for the = 10 results. We suspect two 
effects playing against a fast convergence: -i) the use of frozen spectroscopic orbitals 
fixed to the Hartree-Fock solution, -ii) the independent optimization of the PCF orbital 
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Figure 7. Convergence of the hyperfine contact parameter with respect to the largest 
principal quantum number of the active sets for the ground state of lithium. The 
reference value (dotted line) correspond to the results obtained by Yan et al. [A2\ . 



Table 2. Energies, ^sms E^nd Of^ont increasing maximum principal quantum 

number for the ground state of Li. 



CAS-CP-DPCFI {Is, 2s} 





Energy (a.u.) 


5'sms 


•^cont 


HF 


-7.432726927 


0.000000000 


2.0932317 


4 


-7.476750919 


0.304445112 


2.8624693 


5 


-7.477417364 


0.303171470 


2.8984561 


6 


-7.477689394 


0.302588400 


2.8996669 


7 


-7.477824739 


0.302322528 


2.9007289 


8 


-7.477900177 


0.302162573 


2.9019278 


9 


-7.477946240 


0.302062836 


2.9025948 


10 


-7.477975286 


0.302006041 


2.9033303 


Yan et al. gS 


-7.47806032310(31) 


0.301842809(15) 


2.905922(50) 



sets that forbids the couphng between the PCFs subspaces to capture the indirect effects 
in the orbital optimization. The unconstrained solution could be obtained by solving a 
general non-orthogonal MCHF problem mixing the different PCFs for the orbital SCF 
optimization, and substituting, at the end of each orbital self-consistent-field step, the 
conventional CI by a DPCFI (ie. non-orthogonal CI) calculation to get the desired 
deconstrained eigenvector. 

7. Partial deconstraining schemes 

It is worthwhile to investigate how indirect effects are distributed over CSFs. To an- 
swer this question, we focus our analysis on the mixing coefficient constraint itself by 
eliminating the impact of the orbital optimization. We are therefore diagonalizing the 
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Figure 8. Convergence of the hyperfine parameters with respect to the largest 
principal quantum number of the active sets for the first excited state, ls^2p 2P° 
of lithium. The reference value (dotted line) corresponds to the result obtained by 
Yan et al. l42l . 



Hamiltonian operator within each PCF space, adopting the same orthonormal orbital 
common set for all PCFs. At each step we promote one selected CSF, included in 
one PCF, at the same level than the reference CSF and then we solve the associated 
eigenvalue problem for finding the total energy and the two other properties. In this 
way, we progressively remove all constraints in the expansion coefficients, going from 
a low-dimension PCFI to a larger CSF-CI calculation. Our selection rule for choosing 
at each step the promoted CSF is somewhat arbitrary since we adopt the order of the 
configuration list produced by Isgen program |46]. In the present case, the hierarchy 
follows i) S from the valence, ii) S from the core, iii) D from core- valence and iv) D from 
the core. Figure |9] illustrates the evolution, for each property, of the relative difference 
between the value calculated for a given matrix size, associated to a partially decon- 
strained many-electron wave function, and the corresponding CSF-CI result. Figure [9] 
shows not only that the operators are affected differently by the constraint, but also that 
the way in which the constraint effect is distributed over the CSFs directly depends on 
the selected operator itself. One may note that, accordingly to the variational principle. 



A PCFI approach for describing electron correlation in atoms 



23 



Table 3. Energies, S'sms and hyperfine parameters obtained using the PCFI method 
for the ^-P° state. 



CAS-CP-DPCFI {Is, 2p} 



'^max 


Energy (a.u.) 


S'sms 


"orb 


HF 


-7.365069658 


-0.041898309 


0.0585716 


4 


-7.408715289 


0.249999650 


0.0611674 


5 


-7.409489337 


0.248269629 


0.0626329 


6 


-7.409782422 


0.247618411 


0.0627899 


7 


-7.409921067 


0.247249783 


0.0628740 


8 


-7.409996500 


0.247032302 


0.0629599 


9 


-7.410041507 


0.246907864 


0.0629917 


10 


-7.410071044 


0.246830104 


0.0630317 


Van pf nl 1411 

1 Clli tit. 1^11 




n 0Af^'7'^'^f^^ ^7^ ^ 






^dip 


<^cont 


^quad 


HF 


-0.0117143 


0.0000000 


-0.0234287 


4 


-0.0130095 


-0.2078035 


-0.0220013 


5 


-0.0133191 


-0.2122187 


-0.0225979 


6 


-0.0133551 


-0.2122960 


-0.0226542 


7 


-0.0133691 


-0.2128985 


-0.0226978 


8 


-0.0133940 


-0.2130121 


-0.0227275 


9 


-0.0133960 


-0.2133731 


-0.0227428 


10 


-0.0134106 


-0.2134760 


-0.0227578 


Yan et al. 42; 




-0.214860(50) 





Table 4. Line strength and oscillator strength for the — ^ S resonance line. 



CAS-CP-DPCFI {Is, 2/} 



^max 


AE (cm-i) 


Si 


Sy 




gfi 


gfv 


4 




14930.92 


33.4862 


32 


7741 


1.51872 


1.48642 


5 




14907.31 


33.1284 


32 


9457 


1.50011 


1.49184 


6 




14902.69 


33.0804 


32 


9798 


1.49748 


1.49292 


7 




14901.97 


33.0639 


32 


9835 


1.49666 


1.49302 


8 




14901.97 


33.0536 


32 


9761 


1.49619 


1.49268 


9 




14902.20 


33.0465 


32 


9707 


1.49589 


1.49246 


10 




14902.09 


33.0414 


32 


9692 


1.49565 


1.49238 


CAS (n= 


10) 


14902.23 


33.0038 


32 


9973 


1.49396 


1.49367 


CAS (n= 


10) {Is, 2/} 


14902.48 


33.0075 


32 


9885 


1.49415 


1.49329 


Yan et a 


• m 


14903.16176(29) 


33.00066933 


33 


00081733 


1.4939139 


1.4939206 



the total energy is monotonically decreasing until it reaches the CI value, but the evolu- 
tion is much more erratic for the two other properties. The presence of many plateaus 
in the graphs of figure [9] suggests that some efficient partial deconstraint scheme might 
be found. For the contact term, the latter is accidentally revealed, as illustrated by the 
third graph of figure [HI due to the fact that the most important contributions arise from 
the the single excitations from the core that appear first in the generated list. For the 
S'sms parameter, the jump occurring around 400 corresponds to the presence of 2s2p^ "^S 
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Figure 9. Progressive deconstraint of the many-electron wave function. The 
three curves illustrate the relative difference between the CI calculation and the 
partially deconstrained problem for respectively the total energy, the specific mass 
shift parameter and the hyperfine contact paremeter for the ground state ls^2s '^S of 
neutral lithium. 



in the expansion of Ais_isis. This excitation is introduced quite late in the sequence 
since it is generated by Isgen in the fourth category (double excitations from the core). 

By selecting and deconstraining the CSFs which constitute the dominant 
contributors to the constraint effect, it should be possible to optimize the CSF 
space partition that produces the more efficient PCF basis. We will show that full 
deconstraining is not always necessary and that some deconstraining schemes might be 
better than others. The key point in this analysis is to find in the CP space, for a 
given property, the best candidates for a possible promotion in the CSF transfer ( IT71) . 
As illustrated by figure [H it is hard to predict the behavior of a given property with 
respect to a particular scheme of deconstraint, except for the total energy, thanks 
to the Hylleraas-Undheim theorem |17]. It is here interesting to remember what we 
learn from the time- independent perturbation theory. Defining the \n^'^^) states as the 
eigenstates of the zeroth-order Hamiltonian H^^\ and introducing a perturbation V 
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such as H = + AV^, the eigenfunction \n) of H can be expressed as the following 
A-expansion 

l"> = l"'°'> + ^El*™>TJi5^ (29) 

kj^n ^k 

\k^''^)VkiVin ^ |A:(°))V;„\4n \ 

|y|j(o)^ matrix elements are calculated using the zero-order 
basis [H] . From this expression, it is easy to derive the A-expansion of the expectation 
value of any operator X 

{n\X\n) = (n(°)|X|n(°)) + A f V , + ^ 



'I] 



(0) 



, 2 ( \^ Xkk'VknVk'n 

\k^nk'^n\^ri. r,}. ){J^n -C/^/J 

XnkyklVln XnkynnVkn 

f^(o) _ r(o)vr(o) _ f^(0) _^(0))2 

fc^n i^n l-^'^ -^A: J l"^" "^i J /cT^n l"^" "^fc J 
\ ^ \ ^ XknVklVln \ ^ XknVnnVkn \ 

The first term is the reference value. The 0{X) terms are called the "crossed second- 
order" (cso) corrections to the zeroth-order expectation value {X)n since they involve 
both the Hamiltonian and the operator X in the coupling matrix elements. These 
cso contributions have been used successfully in the parametric method to analyze 
isotope shifts |19], hyperfine structures [50] and field shifts [51]. From this expression, 
one can realize that it would be interesting for an efficient ab initio approach, to 
partially deconstrain the PCFs according to their contribution to the expectation value. 
The CSFs that are the best candidates for a transfer promotion (ITTll are the ones 
that are coupled to the reference by both the Hamiltonian and the relevant operator 
{Xnk 7^ 0; Vnk 7^ 0) but these are not the only ones. For the SMS parameter of the 
lithium ground state, the HF value is strictly zero. The non zero value is made of 
correlation components only 

Ssms = 5^QC,(<l>,|Xs^s|$i) = 0.301 450 504 (31) 

ij 

and a close ranking analysis of the SMS matrix shows that the first ten contributions 
to ( 1311) reported in Table [5] account for 95% of the total contributions. 

In a perturbation approach, the contributions appearing in lines (1, 3, 6, 7, 8) 
of table |5] would be cso contributions of 0{X). The other five higher-order corrections 
are rather important too and would result from contributions of the type (Xkk'VknVk'n)- 
If we use, for the progressive PCF deconstraint, the list of the CSFs sorted according to 
their contribution to the Ssms parameter we get a very efficient deconstraining scheme 
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Table 5. The first ten more important contributions to ([31]) for the hthimn gromid 
state S'sms parameter. 
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Figure 10. Li ls^2s ^5*. Partial deconstraint according to the SMS contributions 
sorting. 



for this property. Figure [TO] illustrates indeed that the transfer of the ~ 140 first CSFs 
produces a rather good result, to put in contrast with Figure [9] 

For the total energy, the most natural deconstraining scheme is based on a the 
weight criteria. The importance of a CSF is indeed mainly determined by the magnitude 
of its mixing coefficient that can be estimated from the first-order correction to the wave 
function (the 0{\) term in f l29]) ). i.e. 

~ p(o) (0) (32) 

To illustrate the gain of partial deconstraint, we focus on Be ls^2s2p considered 
in section [6^ Starting from the constrained representation corresponding to the PCFI 
solution, we deconstrain, in each of the three MR-PCF wave function f ll2p . all CSFs 
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having a mixing coefficient larger than a given threshold, ie. |c^| > e. Repeating this 
operation for a decreasing threshold e, we progressively move from the constrained PCFI 
solution (e = 1) to the deconstrained DPCFI one (e = 0). Using e = 1. 10""^, we reduce 
the computer time by a factor of 10 — 20 relatively to the DPCFI calculations. For 
the n = 10 active set, the size of the eigenvalue problem is respectively 9 (PCFI), 788 
(partially deconstrained) and 21 632 (DPCI) wave functions. The values for the total 
energy, the S^ms and hyperfine parameters are compared with the HF, MR, [n = 10) 
SD-MR-MCHF, PCFI and DPCFI results in Table El The (n = 10) PCFI result is 



Table 6. Comparison between constrained PCFI, partially deconstrained (e = 1. 10 
and DPCFI results for Be ls^2s2p ^P°. 
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MR 
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deconstrained SD-MR-PCFI 


9 
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10 
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SD-MR-DPCFI 


10 
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already better than the {n = 10) SD-MR-MCHF value for the total energy. Taken 
the DPCFI results as the reference for the other parameters, decontraction becomes 
important as discussed in section H] but partial deconstraint with e = 1. 10~^ allows to 
produce quite accurate results, recapturing 96 % (S'sms)? 93 % (ctorft) and 52 % {hquad) 
of the way to the DPCFI values. The fact that only 52 % of hquad is recaptured once 
again illustrates the weak coupling of this quantity to the energy. 

8. Conclusions and perspectives 

We extend the multi-configuration Hartree-Fock method by partitioning the correlation 
function space into several subspaces to get a better description of the dynamical 
correlation. For each of the subspaces, we optimize a partition correlation function to 
correct the multireference wave function. The atomic state is described by the eigenstate 
of a low-dimensional PCFI eigenvalue problem, regrouping the different PCFs with the 
multireference function. By relaxing the orthonormality restriction on the orbital basis 
between the different subspaces, one obtains several mutually non-orthogonal orbital 
basis sets that are better adapted to the short range nature of the dynamical correlation. 

Our original work [20], targeting the ground state of beryllium, leads to some 
unexpected complication when the PCFI eigenvector is used for the calculation of 
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spectroscopic properties other than the total energy. We identified the source of the 
problem as the loss of variational freedom in solving the interaction eigenvalue problem 
that determines the PCF weights. We evaluated the differences in terms of convergence 
between properties calculated, on the one hand with a single orthonormal basis spanning 
an MCHF wave function in the traditional approach, and on the other hand with several 
tailored orbital basis sets optimized for different PCFs and underlying a (D)PCFI wave 
function. In addition to energy we consider specific mass shift, hyperfine structure 
parameters and transition data. We demonstrate that by (partially) deconstraining the 
mixing coefficients within each PCF subspace, one converges to the correct limits and 
keeps the tremendous advantage of improved convergence rates that comes from the 
use of several orbital sets. Reducing ultimately each PCF to a single CSF with its own 
orbital basis leads to a non-orthogonal configuration interaction approach. 

We found that partial deconstraining schemes are very attractive. For the energy 
indeed, the PCFI approach is shown to be highly efficient by producing accurate 
eigenvalues. Partial deconstraining schemes based on the CSF weight criteria appear 
efficient and natural for the Hamiltonian. For other properties, we demonstrated that a 
deconstraint scheme based on the ranking of the contributions to the relevant observable 
becomes much more efficient, as expected from a perturbation analysis. The PCFI 
variational method, assisted by such a pertiubation analysis to define the CSFs that 
should be targeted in the deconstraining strategy, is extremely promising in view of 
the present results. The interest of the method is the fast convergence of the desired 
observable combined with a large reduction factor in the size of the interaction problem 
(and the expected CPU time gain factor that goes with it), added to the possibility of 
parallelizing the PCF-optimization processes before regrouping the various subspaces 
with the reference function to build the atomic state. 

The (D)PCFI method is also highly flexible. Taking advantage of the possibility 
of tailoring the correlation subspace partitions, we built a specific partition correlation 
function for capturing core-polarization effects. The addition of this PCF tremendously 
improves the convergence of all the hyperfine parameters. This approach provides values 
of the total energy, specific mass shift and hyperfine parameters with an acceptable 
accuracy for rather small (n = 6) orbital active sets. These calculations demonstrate 
that it becomes easy to gain freedom for capturing effects weakly connected to energy, 
improving the convergence of atomic properties other than the energy. 

For the three- and four-electron systems considered in the present work, the 
(D)PCFI method might not appear as really profitable since rather complete MCHF 
calculations with a single orthonomormal orbital basis can be still performed. For 
larger systems however, as explained in the introduction, the scaling wall becomes the 
real barrier and there is no hope to capture the dynamic correlation between electrons 
in all the different shells. The breaking down of the very large calculations based on a 
common orbital set into a scries of smaller parallel calculations will then become really 
relevant and opens new perspectives. With this respect, five- and six-electron systems 
that are beyond the current hmits of the Hylleraas method are the first targets in our 
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research list. 

The proposed (D)PCFI method can be apphed to any multiconfiguration approach. 
We will soon start to implement it in the fully relativistic codes [52] but the present 
calculations in the non-relativistic framework [53] remain the ideal benchmarks for 
developing new computational strategies. Another natural line of research is to 
investigate how the remaining constraints on radial orbitals could be fully relaxed. 
This long term project could lead to a more general non-orthogonal MCHF approach 
corresponding to the unconstrained SCF-DPCFI approach. 

In quantum chemistry, a method for a one-component relativistic treatment of 
molecular systems using a multiconfigurational approach (CASSCF) with dynamical 
correlation treated with second order perturbation theory (CASPT2) has been proposed 
by Malmqvist et al [55]. In their work, the Hamiltonian matrix is obtained by 
an extension of the restricted active space state interaction (RASSI) method. They 
also use non-orthogonal orbitals and open some interesting perspectives that we might 
investigate for future work, combining RASSI and the present PCFI approach. 
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Appendix 1: The PCF building rules for a unit PCFs overlap matrix. 

Let us assume that the MR function is built on the reference orthonormal subset {(pk} '■ 

and consider two different PCFs, Ai ct A2, built on the MR orbitals and on independent 
one-electron orbital sets {9k} and {Cfe}, 

• |Ai) using an orthonormal orbital basis {0i, 0„ ; ^1, ^mi}, 

• IA2) using anoi/ier orthonormal orbital basis {0i, ...,0ri ; Ci, ■■■■iCm2\ , 

both built from excitations of the ( mult i-) reference function \[''^^. In our approach, the 
reference orbitals subsets {0} of both PCFs Ai and A2 arc necessarily identical, being 
taken from the optimization of the reference function (\1''^''^). 

In second quantization, one can write the orbital orthonormality within the three 
subsets {0}, {9} and {C} as follows 

{(t>a, (t>h] = {(/'L </>!} = and (/)J} ^ 5a,b 

{9a, 9b} = {9l 91} = and {9^, ^J} = 5a,b 
{Ca, Cb} = {Cl Cl} = and {Ca, Cj} = Sa,b . 
The orbital orthogonality being imposed within each PCF, we also have 

{<Pa,9l} = {^aXl} = 0- 

On the other hand, the orthogonality between {9} and {(} orbital subsets is lost : 
{da, Cb} = {01 Cl} = With {9a, cl} = . 

An arbitrary PCF can be written as a linear combination of CSFs and consequently, 
of Slater derminants (Sdets) . We will show that the overlap between any pair of Sdets 
appearing in the overlap matrix element involving two different PCFs is zero. A corollary 
is that the overlap between PCFs is inevitably zero, ie. (Aj|Aj) = 5ij. 
Let 

• {0\(f)i...(f)n---9i...9n , an arbitrary Sdet belonging to PCF (Ai|, 

• Cl---Ci---0n---0i|O) , another Sdet belonging to the other PCF IA2). 
Let us assume that the PCFs are built according the following rules: 

(i) the excitation prototypes appearing in different PCFs have necessarily different 
occupation numbers for the reference orbitals {(pk}, 

(ii) there is no CSF redundancy (Aj fl Aj = 0), 

(iii) the orthogonality constraints between the variational orbitals, {9k} in Ai, and {(k} 
in A2, and the (frozen) reference ones {(/)k} are imposed as usual, through the 
Lagrange multipliers of the MCHF equations. 
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By respecting these ^''building rules" for PCFs constructed from excitations of or- 
bitals {(j)k} spanning \1>^^, we are assured that Sdets belonging to different PCFs have 
different occupation numbers for the reference orbitals {4>k}- This imphes that at least 
one of the creation operators associated with a reference spin-orbital appears only in 
one of the two Sdets, and not in the other one. This creation operator is then free to 
"sail" in the creation/annihilation operator sea, with a possible phase factor as the only 
price to pay. Moving it to the left/right hand-side until it acts on the bra/ket vacuum 
state, we ultimately get a null result (qed) 

As an illustration, if the 02 spin-orbital does not appear in the bra, we have : 

(O|0i03...0n...^l...^nCl...Cl...0l-0^0l|O) 

= (-l)^^(O|050i03...0n...^^l...^^nCl-C|-0l-0l|O) = 

where A''^ is the required number of permutations. 

We have demonstrated that the PCFs overlap matrix is the unit matrix if we respect 
the above PCF building rules. There are however situations for which violating these 
rules might be interesting from the variational point of view. In these cases, the PCF 
overlap matrix is not the unity matrix anymore, forcing us to treat the PCFI interaction 
problem as the generalized eigenvalue problem (fT^ . 

The non-diagonal PCF overlap matrix element (A1IA2) will indeed differ from zero if 

• there exist two CSFs, belonging respectively to Ai and A2, that have identical 
occupation numbers for all the reference orbitals, 

• the reference orbitals cannot be strictly qualified as "reference" ones, when another 
subset is used for one of the PCFs ie. {0^} 7^ {0fc} • For example, if the {0'^} 
orbitals of Ai do not have the same radial functions than the {(f)k} of and A2, 
although having the same labels, some commutation rules are lost : 

• non-orthogonalities are introduced within a given PCF, for instance within A2: 
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Appendix 2: The weight matrix formahsm 



In the present appendix, we will show how the whole sub-matrix block of the 
matrix f|T5|) can be evaluated by performing a single biorthonormal transformation 
treating simultaneouly the counter-transformation of all the elements constituting the 
block-basis. For this purpose, we introduce the weight matrix D*, corresponding to the 
deconstrained PCF {^\, . . . , Aj }. Each column of the latter is composed of the CSF 
expansion coefficients for each of the {hi + 1) elements 
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The number of lines of this weight matrix is given by the size of the original PCF 
{dim{Ai)) to which we add dim{Cl]Di) CSFs necessary for satisfying the CUD condition. 
In the deconstrained limit case, the weight matrix is a rectangular unit matrix of size 
([dim(Ai) + dim{C\JDi)] x dim(Ki)) 
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Adopting the weight matrix, the expression of a sub-matrix block {{hi -|- 1) x {hj + 1)) 
involving two different deconstrained PCFs can be obtained by a simple matrix 
multiplication 

. . . , ^\X}\Om^ • • • , = {BjOB^ (.3) 

involving the counter-transformed weight matrix 
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a similar one for D-' , and the matrix representation of the selected operator O f l2ip . The 
weight matrix formalism introduces the possibility of doing the simultaneous counter- 
transformation, associated to the biorthonormal transformation, of the expansion 
coefficients of each element of the deconstrained PCF. As each column of the original 
weight matrix ([T]) represents the expansion coefficients of each element, i.e. . . . , 
and A^, of a deconstrained PCF in the original CSF space, each column of the counter- 
transform weight matrix ([3]) gives the representation of these elements in the new CSF 
space expressed in the biorthonormal orbital basis. Another advantage of the weight 
matrix formalism is that we are able to treat, without any distinction, an arbitrary 
degree of deconstraint. The price to pay is a simple matrix dimension adaptation. By 
allowing high degrees of deconstraint, we quickly get large weight matrices and the use 
of their sparse representation is welcome. 

For simplifying the counter-transformation process associated to the biorthonormal 
transformation, the Iscud program originally developed for closing a CSF expansion 
under de-excitation [20] was adapted to provide the sparse structure of the counter- 
transformed weight matrix D* before performing the biorthonormalisation by itself. This 
is made possible by the use of the Warshall's algorithm [M] applied to the adjacency 
matrix mapping the single de-excitation operator for all /-values within the CSF basis. 
The obtained transitive closure matrix of the single de-excitation operator corresponds 
to the mask of the D* matrix that determines the position of the non-zero elements after 
biorthonormalisation. 
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